home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / MFSD.FOR < prev    next >
Text File  |  1985-11-29  |  4KB  |  110 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE MFSD
  5. C
  6. C        PURPOSE
  7. C           FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
  8. C
  9. C        USAGE
  10. C           CALL MFSD(A,N,EPS,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           A      - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC
  14. C                    POSITIVE DEFINITE N BY N COEFFICIENT MATRIX.
  15. C                    ON RETURN A CONTAINS THE RESULTANT UPPER
  16. C                    TRIANGULAR MATRIX.
  17. C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
  18. C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
  19. C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  20. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  21. C                    IER=0  - NO ERROR
  22. C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
  23. C                             TER N OR BECAUSE SOME RADICAND IS NON-
  24. C                             POSITIVE (MATRIX A IS NOT POSITIVE
  25. C                             DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-
  26. C                             FICANCE)
  27. C                    IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFI-
  28. C                             CANCE. THE RADICAND FORMED AT FACTORIZA-
  29. C                             TION STEP K+1 WAS STILL POSITIVE BUT NO
  30. C                             LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).
  31. C
  32. C        REMARKS
  33. C           THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE
  34. C           STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.
  35. C           IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-
  36. C           LAR MATRIX IS STORED COLUMNWISE TOO.
  37. C           THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL
  38. C           CALCULATED RADICANDS ARE POSITIVE.
  39. C           THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE
  40. C           SQUARE-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX.
  41. C
  42. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  43. C           NONE
  44. C
  45. C        METHOD
  46. C           SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY.
  47. C           THE GIVEN MATRIX IS REPRESENTED AS PRODUCT OF TWO TRIANGULAR
  48. C           MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF
  49. C           THE RETURNED RIGHT HAND FACTOR.
  50. C
  51. C     ..................................................................
  52. C
  53.       SUBROUTINE MFSD(A,N,EPS,IER)
  54. C
  55. C
  56.       DIMENSION A(1)
  57.       DOUBLE PRECISION DPIV,DSUM
  58. C
  59. C        TEST ON WRONG INPUT PARAMETER N
  60.       IF(N-1) 12,1,1
  61.     1 IER=0
  62. C
  63. C        INITIALIZE DIAGONAL-LOOP
  64.       KPIV=0
  65.       DO 11 K=1,N
  66.       KPIV=KPIV+K
  67.       IND=KPIV
  68.       LEND=K-1
  69. C
  70. C        CALCULATE TOLERANCE
  71.       TOL=ABS(EPS*A(KPIV))
  72. C
  73. C        START FACTORIZATION-LOOP OVER K-TH ROW
  74.       DO 11 I=K,N
  75.       DSUM=0.D0
  76.       IF(LEND) 2,4,2
  77. C
  78. C        START INNER LOOP
  79.     2 DO 3 L=1,LEND
  80.       LANF=KPIV-L
  81.       LIND=IND-L
  82.     3 DSUM=DSUM+DBLE(A(LANF)*A(LIND))
  83. C        END OF INNER LOOP
  84. C
  85. C        TRANSFORM ELEMENT A(IND)
  86.     4 DSUM=DBLE(A(IND))-DSUM
  87.       IF(I-K) 10,5,10
  88. C
  89. C        TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE
  90.     5 IF(SNGL(DSUM)-TOL) 6,6,9
  91.     6 IF(DSUM) 12,12,7
  92.     7 IF(IER) 8,8,9
  93.     8 IER=K-1
  94. C
  95. C        COMPUTE PIVOT ELEMENT
  96.     9 DPIV=DSQRT(DSUM)
  97.       A(KPIV)=DPIV
  98.       DPIV=1.D0/DPIV
  99.       GO TO 11
  100. C
  101. C        CALCULATE TERMS IN ROW
  102.    10 A(IND)=DSUM*DPIV
  103.    11 IND=IND+I
  104. C
  105. C        END OF DIAGONAL-LOOP
  106.       RETURN
  107.    12 IER=-1
  108.       RETURN
  109.       END
  110.